

set.seed(500)

make_simulation_did <- function(tau, sd, alpha = 0.05, sims = 250) {
  
  possible_ns <- seq(from = 100, to = 2000, by = 100)
  powers <- rep(NA, length(possible_ns))
  
  for (j in seq_along(possible_ns)) {
    
    N <- possible_ns[j]
    significant_experiments <- rep(NA, sims)
    
    for (i in 1:sims) {
      
      treated_n <- round(0.2 * N)
      control_n <- N - treated_n
      
      # Each unit has pre and post observations (repeated cross-section structure)
      # So total rows = 2 * N
      total_n <- 2 * N
      
      group <- c(rep(0, 2 * control_n), rep(1, 2 * treated_n))  # control = 0, treated = 1
      time <- rep(c(0, 1), times = N)  # pre = 0, post = 1
      
      interaction <- group * time
      
      error <- rnorm(total_n, mean = 0, sd = sd)
      Y <- 2.8 + tau * interaction + error
      
      dat <- data.frame(Y = Y,
                        pre_post = time,
                        validated_election_outcome = group)
      
      fit <- try(lm_robust(Y ~ pre_post * validated_election_outcome, data = dat), silent = TRUE)
      
      if (inherits(fit, "try-error")) next
      
      p_value <- summary(fit)$coefficients["pre_post:validated_election_outcome", "Pr(>|t|)"]
      significant_experiments[i] <- (p_value <= alpha)
    }
    
    powers[j] <- mean(significant_experiments, na.rm = TRUE)
  }
  
  return(data.frame(N = possible_ns, Power = powers, tau = tau, sd = sd))
}
results_1_35 <- make_simulation_did(tau = 0.05, sd = 0.20)
results_1_45 <- make_simulation_did(tau = 0.05, sd = 0.30)
results_1_55 <- make_simulation_did(tau = 0.05, sd = 0.40)
results_1_65 <- make_simulation_did(tau = 0.05, sd = 0.50)

results_2_35 <- make_simulation_did(tau = 0.1, sd = 0.20)
results_2_45 <- make_simulation_did(tau = 0.1, sd = 0.30)
results_2_55 <- make_simulation_did(tau = 0.1, sd = 0.40)
results_2_65 <- make_simulation_did(tau = 0.1, sd = 0.50)

results_3_35 <- make_simulation_did(tau = 0.15, sd = 0.20)
results_3_45 <- make_simulation_did(tau = 0.15, sd = 0.30)
results_3_55 <- make_simulation_did(tau = 0.15, sd = 0.40)
results_3_65 <- make_simulation_did(tau = 0.15, sd = 0.50)

results_4_35 <- make_simulation_did(tau = 0.2, sd = 0.20)
results_4_45 <- make_simulation_did(tau = 0.2, sd = 0.30)
results_4_55 <- make_simulation_did(tau = 0.2, sd = 0.40)
results_4_65 <- make_simulation_did(tau = 0.2, sd = 0.50)


results_all <- bind_rows(
  results_1_35 = results_1_35,
  results_1_45 = results_1_45,
  results_1_55 = results_1_55,
  results_1_65 = results_1_65,
  
  results_2_35 = results_2_35,
  results_2_45 = results_2_45,
  results_2_55 = results_2_55,
  results_2_65 = results_2_65,
  
  results_3_35 = results_3_35,
  results_3_45 = results_3_45,
  results_3_55 = results_3_55,
  results_3_65 = results_3_65,
  
  results_4_35 = results_4_35,
  results_4_45 = results_4_45,
  results_4_55 = results_4_55,
  results_4_65 = results_4_65,
  
  .id = "sim_id"
)

p <-
  results_all %>%
  mutate(sd = paste0("SD(Y(0)) = ", sd)) %>%
  ggplot(aes(x = N, y = Power*100, color = factor(tau))) +
  geom_point(size = 1) +
  geom_line() +
  #geom_vline(xintercept = 3000, color = "red", linetype = "dashed") +
  geom_vline(xintercept = 800, color = "pink", linetype = "dashed") +
  geom_hline(yintercept = 80, color = "pink", linetype = "dashed") +
  xlim(0, 2000) +
  #geom_hline(yintercept = 80, color = "red", linetype = "dashed") +
  theme_bw() +
  facet_wrap(sd~., ncol =2) +
  ylab("Probability of Detecting Effect (%)") +
  xlab("Number of Observations") +
  labs(color=expression(beta[3])) +
  scale_color_grey() +
  theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = .2)) + 
  theme(strip.background =element_rect(fill="white")) +
  theme(panel.border = element_blank(),
        axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.line.y = element_line(size = 0.5, linetype = "solid", colour = "black"))

ggsave(filename = "./outputs/figures/figure_a2.pdf", plot = p, width = 8, height = 5)
ggsave(filename = "./outputs/figures/figure_a2.eps", device = "eps", plot = p, width = 8, height = 5)


estimation_data %>% filter(A == 1) %>% pull(avg_score_democ_support) %>% sd(., na.rm = T)
